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Abstract. Unraveling the mechanism for core-collapse supernova explosions is an outstanding 
computational challenge and the problem remains essentially unsolved despite more than four 
decades of effort. However, much progress in realistic modeling has occurred recently through the 
availability of multi-terafiop machines and the increasing sophistication of supernova codes. These 
improvements have led to some key insights which may clarify the picture in the not too distant 
future. Here we briefly review the current status of the three explosion mechanisms (acoustic, MHD, 
and neutrino heating) that are currently under active investigation, concentrating on the neutrino 
heating mechanism as the one most likely responsible for producing explosions from progenitors 
in the mass range ~ 10 to <~ 25 M Q . We then briefly describe the CHIMERA code, a supernova 
code we have developed to simulate core-collapse supernovae in 1, 2, and 3 spatial dimensions. We 
finally describe the results of an ongoing suite of 2D simulations initiated from a 12, 15, 20, and 25 
M progenitor. These have all exhibited explosions and are currently in the expanding phase with 
the shock at between 5,000 and 10,000 km. We finally very briefly describe an ongoing simulation 
in 3 spatial dimensions initiated from the 15 M progenitor. 
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INTRODUCTION 

Much progress has been made unraveling the core collapse supernova mechanism in the 
past ten or so years, but it still remains an unsolved problem. It is apparent from observa- 
tions that core-collapse supernovae exhibit large-scale anisotropics. Spectropolarimetry, 
the large average pulsar velocities, and the morphology of highly resolved images of SN 
1987 A all suggest that anisotropy not only develops, but likely develops very early on 
in the explosion [e.g., see 4, 5, 6, for reviews and references]. This observed asymme- 
try is now believed to be a manifestation of the underlying essential multidimensional 
nature of the supernova mechanism. Realistic numerical core-collapse supernova mod- 
eling thus requires multidimensional techniques and correspondingly massive amounts 
of computer resources to capture this important macro-physics. Additionally, during and 
following core-collapse a plethora of microphysics comes into play and must be com- 
puted accurately as the supernova mechanism appears marginal. Inaccuracies in any part 
of a numerical simulation can prejudice the outcome. It is therefore not surprising that 



the supernova mechanism is taking a long time to unravel. 



MECHANISMS 

Three explosion mechanisms have been the focus of current research: (1) the acoustic 
mechanism, (2) the MHD mechanism, (3) and the neutrino heating mechanism. 



Acoustic Mechanism 

The acoustic mechanism was discovered by Burrows et al. [7] in the simulation 
of an 11 M & progenitor. They found that long after shock stagnation ( >,0.6 s post- 
bounce) turbulence and anisotropic accretion on the proto-neutron star excites and 
maintains vigorous g-mode oscillations which radiate intense sound waves, the energy 
coming from the gravitational binding energy of the accreted gas. As these sound waves 
propagate outward through the negative density gradient into the surroundings, they 
steepen into shocks and their energy and momentum are efficiently absorbed, powering 
up the supernova explosion. Thus, the proto-neutron star acts like a transducer converting 
the gravitational energy of infall into acoustic energy which propagates out and deposits 
energy in the surroundings, in analogy with the neutrino transport mechanism. They 
subsequently found [8] that the acoustic mechanism is able to explode a variety of 
progenitors at late times (,>0.6 s). The physical reality of this mechanism is being 
debated as it has not been observed in simulations by other investigators when these 
simulations have been pushed to late times, although their numerical techniques, though 
different, are capable of capturing this phenomenon. A further note of caution is cast by 
a recent study by Weinberg and Quataert [9] that finds that the damping of the primary 
mode by the parametric instability causes the primary 1 = 1 g-mode to saturate at an 
energy two orders of magnitude lower than that required to power a supernova. 



MHD Mechanism 

Magnetic fields threading a progenitor are frozen in the gas on all relevant core- 
collapse time scales. On core collapse these magnetic fields will be amplified both 
by flux conservation during matter compression and by being wound up toroidally by 
the differential rotation of the core. Simulations with increasing sophistication have 
shown that if the iron core before collapse is threaded by very strong magnetic fields 
(B > 10 12 gauss), then this in combination with rapid rotation can produce jet-like 
explosions magnetically on a prompt time scale [10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. 
Furthermore, it has been recognized that initially weak magnetic fields can be amplified 
to equipartition values exponentially by the magnetorotational instability [20, 21, 22]. 
Notwithstanding all this, it must be appreciated that the maximum magnetic energy that 
can be achieved in a differentially rotating core is the free energy, tf ree , of the differential 
rotation, i.e., the difference between the energy of the differentially rotating core and the 
same core uniformly rotating with the same angular momentum, and 



[22]. Thus rather small initial rotation periods, < 2 ms, for newly formed neutron stars 
are required if enough magnetic energy is to potentially arise to power up the typical 
supernova. These small rotation periods are at variance with the calculated rotational 
periods of the magnetized cores of supernova progenitors [23], and the extrapolated 
periods of newly formed neutron stars (> 10 ms). Both of these constraints are "soft" 
(stellar evolutionary calculations with rotation and magnetic fields are not ab initio, and 
we have not yet observed a newly formed neutron star), but if they hold then the MHD 
mechanism will only be relevant to a subset of core collapse supernovae. However, 
the observations of magnetars, long-duration gamma-ray bursts, and hints of highly 
collimated material in some supernova remnants suggests a subclass of events that are 
magnetically driven. 



Neutrino Heating Mechanism 




FIGURE 1. Entropy and velocity configuration snapshots of the Model 1 lM_Sym_32R. 

The neutrino heating mechanism has a long pedigree extending back to the seminal 
paper of Colgate and White [24] and its more modern incarnation [25]. Following the 
collapse and bounce of the inner core of massive star at the endpoint of its normal 
thermonuclear evolution, the shock launched at core bounce stalls in the outer core, 
losing energy (and therefore post-shock pressure) to nuclear dissociation and electron 
neutrino losses. Within a short time (~ 50 ms) a thermodynamic profile is established 
as illustrated in Figure 1 in which infalling matter encountering the outward flow of 



neutrinos undergoes net heating between the shock and the so-called gain radius, and net 
cooling below, due to the different neutrino heating and cooling radial profiles. Crudely 
speaking, for neutrino heating to be successful in powering an explosion a fluid element 
must be heated sufficiently while it resides in the heating layer to reenergize the shock. 
Sophisticated analyses [26, 27] have shown that the criteria for generating explosions 
in 1-D by neutrino heating depends on both the efficiency of heating between the gain 
radius and the shock and of cooling between the neutrino sphere and the gain radius. 

Energy deposition by neutrinos plays the primary role in the neutrino heating mecha- 
nism, and the rate of energy deposition per nucleon, q, can be written as 

where the first and second terms express the absorption of electron neutrinos (v e 's) and 
antineutrinos (V e 's), respectively. For the v e 's (v e 's), L Ve (Ly e ) is their luminosity, (gy ) 
((e? e )) their mean square energy, and and ( j-) their inverse flux factor, which is a 
measure of their anisotropy. Clearly, an accurate calculation of q requires an accurate 
calculation of both the energy spectrum and the angular distribution of the neutrinos. 



STATUS OF THE NEUTRINO HEATING MECHANISM 

Simulations of core-collapse supernovae in spherical symmetry with considerable real- 
ism have been performed with Boltzmann neutrino transport, state-of-the-art neutrino 
interactions, and with/without general relativity [28, 29, 30, 31]. These have not yielded 
explosions. Something is clearly missing. 

An insight as to what the missing ingredient might be was developed during the 1990's 
and is the essential role played by multidimensional effects. Analyses of immediate post- 
bounce core profiles given by computer simulations had for a long time indicated that 
a variety of fluid instabilities are present [32, 33, 34, 35]. The most important of these 
for the neutrino heating mechanism is the neutrino heating above the neutrino sphere. 
Because neutrinos heat the bottom o f the heating layer most intensely a negative entropy 
gradient builds up which renders the layer convectively unstable. In order for convection 
to grow, however, the fluid must remain in the heating layer for a critical length of time; 
roughly the ratio of the advective timescale to some averaged timescale of convective 

growth timescale must be ~ 3 [36]. If convection can get established in the hearing 
layer, hot gas from the neutrino-heating region will be transported directly to the shock, 
while downflows simultaneously will carry cold, accreted matter to the layer of strongest 
neutrino heating where a part of this gas, being cold, readily absorbs more energy from 
the neutrinos. The loss of energy accompanying the advection of matter through the 
gain radius is thereby reduced and more energy stays in the heating layer. It has been 
shown that convection in the heating layer will lead to explosions where the same model 
with the same neutrino luminosities and RMS energies fail when computed in spherical 
symmetry [37]. 

Another important ingredient missing in spherical symmetry was pointed out by 
Blondin et al. [38] who discovered that the stalled shock is subject to low-mode as- 



pherical oscillations, which they referred to as the standing accretion shock instability 
or 'SASL' The cause of this instability is still being debated Foglizzo et al. [39], Scheck 
et al. [40] arguing for an advective-acoustic cycle and [41, 42] for a purely acoustic cy- 
cle. The development of the SASI leads to an enlargement of the heating layer in one 
region and its diminution in another (Figure 1). Where the heating region is enlarged, 
the efficiency of neutrino heating is enhanced and the onset of convection is made more 
favorable if it has not already commenced. Where the hearing region is constricted, con- 
ditions are favorable for the establishment of down-flows or return-flows for large-scale 
convection. It has further been pointed out and supported by 2D simulations with pa- 
rameterized neutrino sources that the development of the SASI leads to the large asym- 
metries observed for SN 1987 A and other supernovae, and might account for the large 
observed velocities of neutron stars [43, 44]. 

Several groups, the Arizona-Jerusalem collaboration, the Garching group, and the 
FAU-NCSU-Oak Ridge collaboration have modeled core collapse in 2 spatial dimen- 
sions with spectral neutrino transport and the results so far do not exhibit convergence. 
The Arizona- Jeruseleum group have not seen neutrino driven explosions in any of their 
simulations [7, 16]. They apply 2D neutrino diffusion but neglect any energy-bin cou- 
pling and velocity corrections, and employ Newtonian gravity. The Garching group 
[45, 46, 47, 48] apply a "ray-by-ray-plus" treatment of neutrino transport, where trans- 
port is calculated along radial rays, neglecting neutrino shear and nonradial fluxes, but 
including all other lateral effects. Transport is computed by solving the neutrino number, 
momentum, and energy equations closed by a variable Eddington factor computed from 
a simplified Boltzmann equation. The neutrino microphysics is state of the art. Gravity 
is computed using an approximate general relativistic potential for spherical gravity and 
Newtonian higher moments. The Garching group has evolved an 11.2 M and a 15 M 
progenitor and finds explosions commencing after rather long post bounce times, 220 
ms and 620 ms respectively. Their simulations were not carried long enough to ascertain 
the explosion energies directly. Using the CHIMERA code, our group (the FAU-NCSU- 
Oak Ridge collaboration) has carried out simulations for a suite of four progenitors of 
MS masses 12 M , 15 M , 20 M ,and 25 M evolved to core collapse by Woosley 
and Heger [49]. We obtain explosions for all these models and find that they power-up 
earlier compared with the Garching simulations. Our simulations are ongoing with the 
models having currently been evolved to post bounce times of from 400 to 700 ms. The 
CHIMERA code and our results will be described below. 

THE CHIMERA CODE 

The CHIMERA code is designed to simulate core-collapse supernovae in 1, 2, and 
3 spatial dimensions from the onset of collapse to the order of 1 sec post bounce 
given present day state-of-the-art computational resources, such as the Cray XT4. It 
conserves total energy (gravitational, internal, kinetic, and neutrino) to within ± 0.5 B. 
The code currently has three main components: a hydro component, a neutrino transport 
component, and a nuclear reaction network component. In addition there is a Poisson 
solver for the gravitational potential and a sophisticated equation of state. A preliminary 
version of the code was briefly described in Bruenn et al. [50], Mezzacappa et al. [51]. 



The hydrodynamics is evolved via a Lagrangian remap implementation of the Piece- 
wise Parabolic Method (PPM) [1]. A moving radial grid option wherein the radial grid 
follows the average radial motion of the fluid makes it possible for the core infall phase 
to be followed with good resolution. Following bounce, an adaptive mesh redistribution 
algorithm keeps the radial grid between the core center and the shock structured so as to 
maintain approximately constant Ap/p . For 256 radial zones, this ensures that there are 
at least 15 radial zones per decade in density. For 512 radial zones (our higher resolu- 
tion runs) at least 30 radial zones per decade in density are maintained. The equation of 
state (EOS) blend of Lattimer and Swesty [52], highly modified Cooperstein [53], and 
other details are described in Bruenn et al. [50]. The algorithm described in [54], mod- 
ified for greater rubustness, was employed to stabilize shocks oriented along grid lines. 
The gravitational potential is solved by means of the Newtonian gravity spectral Poisson 
solver described in [55], and modified for application to 3D. Gravity was computed as 
described above for the Garching group. 

Neutrino transport is implemented by a "ray-by-ray-plus" approximation [cf. 46] 
whereby the lateral effects of neutrinos such as lateral pressure gradients (in optically 
thick conditions), neutrino advection, and velocity corrections are taken into account, 
but transport is performed only in the radial direction. Transport is computed by means 
of multigroup flux-limited diffusion with a sophisticated flux limiter that has been tuned 
to reproduce results of a general relativistic Boltzmann transport results to within a few 
percent [56]. All 0(v/c) observer corrections have been included. The transport solver 
is fully implicit and solves for four neutrino flavors simultaneously (i.e., v e 's, V e 's, v^'s 
and v T 's (collectively V x 's), and v^'s and v T 's (collectively V x 's)), allowing for neutrino- 
neutrino scattering and pair-exchange, and different V and V opacities. State-of-the-art 
neutrino interactions are included with full energy dependences. 

The nuclear composition in the non-NSE regions of these models is evolved by the 
thermonuclear reaction network of [3]. This is a fully implicit general purpose reaction 
network, however in these models only reactions linking the 14 alpha nuclei from 4 He to 
60 Zn are used. Data for these reactions is drawn from the REACLIB compilations [57]. 
The nucleons have only very small abundances at any time and are included to make the 
NSE-nonNSE transition smoother. The iron-like nucleus is included to conserve charge 
in a freezeout occurring with an electron fraction below 0.5. 

2D SIMULATION RESULTS 

With the many improvements in the CHIMERA code and additions of more realistic 
physics (e.g., GR corrections, state-of-the-art neutrino microphysics) since our first 
results were reported in Bruenn et al. [50] we find that all our models explode as before, 
but now the nuclear energy released when the oxygen layers reach the shock does not 
play the critical role it did before. For a period of time after bounce (65 ms for the 12 
and 15 M models, 105 ms and 130 ms respectively for the 20 and 25 M models) the 
2D and corresponding spherically symmetric simulations track each other very closely. 
(This is apart from a brief episode in the 2D models of lepton driven convection in the 
proto-neutron star, followed by an equally brief episode of neutrino driven convection 
from the negative entropy gradient imprinted by the shock as it weakened.) During this 



FIGURE 2. Two snapshots of the entropy distribution in the 15 M model at times 124 ms 
(left) and 255 ms (right) post-bounce. Note the difference in scale of the two figures, 



period the heating layer becomes convectively unstable but its growth is suppressed in 
the 2D simulations because of the too rapid inflow of material. At the end of this period 
convective mushrooms begin to appear in the heating layer and at almost the same time 
the shock begins to exhibit SASI dipole and quadrupole deformations along the polar 
axis. This expands the heating layer in the polar regions and the enhanced convection 
in these regions quickly develops into a large-scale overturn (Fig. 2, left), with material 
rising and expanding in the polar regions and flowing towards the equatorial region. By 
225 ms post-bounce for the 15 M model (Fig. 2, right), as an example, large lobes of 
high-entropy neutrino heated gas are pushing the shock outwards in the polar regions 
while an equatorial accretion funnel has developed which is channeling lower-entropy 
newly shocked material down to the gain layer. The narrow end of the accretion funnel 
tends to oscillate slowly from one hemisphere to the other, seeming thereby to alternately 
pump up one hemisphere and then the other with energy by the enhanced neutrino 
emission at its base. Simultaneous with this convective activity is a large dipolar SASI 
oscillation of the shock on a more-or-less fixed quadrupolar background. It is likely that 
the lateral sloshing of material accompanying this SASI is what causes the base of the 
accretion funnel to oscillate between hemispheres, i.e. the dog is indeed wagging the tail 
rather than the other way around. 



Eventually a runaway condition is met and the shock begins to accelerate rapidly 
outwards. Measured by the first signs of positive radial post-shock velocities (a different 
criterion than that used by the Garching group) explosions commence at roughly 300 ms 
from bounce for the four models. The 12, 15, and 25 M models exhibit a pronounced 
quadrupolar prolate configuration of the shock as the explosion commences, while the 
20 M the configuration of the shock is more dipolar. The ratio of the major to minor 
axis of the shock ranges from 1.5 (25 M model) to 2.0 (15 M model). The explosion 
energies, computed conservatively as the sum of the kinetic, internal thermal, and kinetic 
energies of outward moving material is shown in Figure 3. 
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FIGURE 3. Explosion energies as a function of post-bounce time. 



3D SIMULATION RESULTS 

We mention finally an ongoing medium resolution simulation in 3 spatial dimensions 
with 304 radial zones, 76 angular, and 152 azimuthal zones, initiated from the 15 M 
progenitor. The simulation has progressed to 139 ms post-bounce and is exhibiting large 
scale convection patterns. We anticipate that this simulation will tell us much about the 
limitations of modeling core collapse in 2D. 
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